(in-ns 'miked.math)
(clojure/refer 'clojure)

(import '(java.util BitSet))

(defn divides 
  "Returns whether n is evenly divisible by d"
  [n d]
  (== 0 (rem n d)))

(defn even
  "Returns whether x is even"
  [x] 
  (divides x 2))

(defn odd [x] 
  (complement (even x)))

(defn get-exp [x base]
  (loop [y 1  exp 0]
    (if (< x y)
      (/ y base)
      (recur (* y base) (inc exp)))))

(defn digits 
  "Returns a list of the digits in the given number"

  ([num base] 
   (defn helper [x y]
     (cons (quot x y)
	   (if (= 1 y)
	     nil
	     (helper (rem x y) (/ y base)))))
   (helper num (get-exp num base)))
  
  ([num]
   (digits num 10))
  
  {:test (fn [] (assert (= '(1 2 3) (digits 123))))})



(defn pow
  "Raise x to yth power"
  {:test (fn [] (assert (= 16 (pow 2 4))))}
    [x y] 
  (if (zero? y) 
    1
    (* x (pow x (dec y)))))

(defn palindrome? [x] 
  (= x (reverse x)))

(defn some-divisor? [x ys] 
  (some (fn [y] (divides x y)) ys))

(defn no-divisor? [x ys] 
  (not-any? (fn [y] (divides x y)) ys))

(defn odd-seq [x] 
  (lazy-cons x (odd-seq (+ x 2))))

(defn square [x] (* x x))

(defn multiples 
  "Sequence of multiples of x"
  {:test (fn [] (assert (= '(3 6 9 12)
			   (take 4 (multiples 3)))))}
  [x]
  (iterate (partial + x) x))

(defn powers
  "Returns a sequence of the powers of the given base" 
  [base] (iterate (fn [x] (* base x))
		  1))

(defn cross-multiply
  {:test (fn [] (assert (= '(3 4 6 8)
			   (cross-multiply '(1 2)
					    '(3 4)))))}
  [xs ys] 
  (mapcat (fn [x] 
	      (map (partial * x) ys))
	  xs))




;;
;; Fib
;;

(defn- fib-helper [a b] (lazy-cons a (fib-helper b (+ a b))))

(def fib
     (rest (rest (fib-helper 0 1))))

;;;
;;; Exponentiation
;;;

(defstruct exp-struct :base :exponent)

(defn exp [base exponent]
  "Returns an exponentiation with base and exponent"
  (struct exp-struct base exponent))


(def base  (accessor exp-struct :base))
(def exponent (accessor exp-struct :exponent))


;;;
;;; Prime cache
;;;


(defn sieve
  "Returns a lazy sieve of prime numbers below limit"

  {:test (fn [] (assert (= '(2 3 5 7 11)
			    (sieve 12))))}
  
  [limit]
  (let [head (ref 2)
	bits (doto (new BitSet (inc limit)) (set 0 limit))
	
	;; When x <= limit, grow the sieve past x and return the new
	;; head, else return nil
	grow (fn []
		 (sync nil
		   (when (> @head 0)
		     (do
		       (doseq y (take-while (partial > limit)
					    (rest (multiples @head)))
			 (. bits (clear y)))
		       (ref-set head (. bits (nextSetBit (inc @head))))
		       (when (> @head 0) @head)))))
	]

    ( (fn produce [] (lazy-cons @head (when (grow)
					(produce)))))))



(defstruct prime-cache-struct :sequence :lookup :limit)


(defn prime-cache [limit]
  (let [s (sieve limit)]

    (struct-map
     prime-cache-struct
     :sequence s
     :lookup (set s) 
     :limit limit)))

(def *prime-cache* (prime-cache 1000000))



(defn prime? 
  "Returns whether x is prime"

  {:test (fn [] (binding [*prime-cache* (prime-cache 10)]
		  (assert (prime? 2))
		  (assert (prime? 3))
		  (assert (prime? 5))
		  (assert (prime? 7))
		  (assert (prime? 11))
		  (assert (prime? 13))))}  

  [x]
  (let [limit (:limit *prime-cache*)]
    (if (< x limit)
      (contains? (:lookup *prime-cache*) x)
      (and (no-divisor? x (:sequence *prime-cache*))
	   (no-divisor? x (take-while (fn [y] (< y x)) 
				      (odd-seq 
				       (if (odd limit)
					 limit (dec limit)))))))))


(defn- next-uncached-prime 
  "Returns the next prime greater than x without referring to the cache"
  [x]
  (if x
    (if (even x)
      (next-uncached-prime (dec x))
      (if (prime? (+ x 2))
	(+ x 2)
	(next-uncached-prime (+ x 2))))
    (next-uncached-prime (:limit *prime-cache*))))


(defn primes 
  "Returns a sequence of prime numbers"
  {:test (fn [] (binding [*prime-cache* (prime-cache 10)]
		  (assert (= '(2 3 5 7 11 13 17 19)
			     (take 8 (primes))))))}
  []
  (concat (:sequence *prime-cache*)
	  (rest (iterate next-uncached-prime nil))))


(defn composite? [x] (not (prime? x)))

;;;
;;; Factoring
;;;

(defn- factor-helper
  ([x primes]
      (let [prime (first primes)]
	   (if (prime? x)  
		(cons x nil)
	      (if (divides x prime)
		  (cons prime (factor-helper (/ x prime) primes))
		(factor-helper x (rest primes)))))))

(defn- bases-to-exps 
  "Given a sorted list of integers, consolidate them into a list of 
   exponentiations."
  {:test (fn [] (assert (= (list (exp 2 3)
				 (exp 3 4))
			   (bases-to-exps '(2 2 2 3 3 3 3)))))}
      
  [factors]
  (let [factor (first factors)  
	f (fn [x] (= x factor))]
    (when factors
      (cons (exp factor (count (take-while f factors)))
	    (bases-to-exps (drop-while f factors))))))

(defn factor
  "Returns the factors of x"
  {:test (fn [] (assert (= (list (exp 2 5)
				 (exp 7 1))
			   (factor 224))))}
    [x] (bases-to-exps (factor-helper x (primes))))


(defn- exp-explode 
  "Returns a list of the base of x raised to all exponents from 1 to the
   exponent of x."
  {:test (fn [] (assert (= '(1024 256 64 16 4 1) (exp-explode (exp 4 5)))))}
  
  [x]
  (reverse (take  (inc (exponent x))
		  (powers (base x)))))

(defn divisors 
  "Returns the divisors of x"
  {:test (fn [] (assert (= '(6 2 3 1) (divisors 6))))}
  [x]
  (reduce cross-multiply (map exp-explode (factor x))))

(defn abundant? [x]
  (< x (apply + (rest (divisors x)))))
